      subroutine Melds3b(filename)
      implicit real*8 (a-h,o-z)
      character*40 fname,filename(20)
      logical filflg
      common /aa/ ia(256000000)
      common /bb/ ib(128000000)
c.... open disp0 file to get the numbers of nodes and degree of freedom
c.... knode .... number of nodes, kdgof .... number of d.o.f.
      open(1,file=filename(1),form='unformatted',status='old')
      read(1) knode,kdgof
      kdeof=kdgof
      close(1)
      inquire(file='time',exist=filflg)
      if (filflg) then
      open(6,file='time',form='unformatted',status='old')
      read(6) tmax,dt,time,it
      close(6)
c      write(*,*) 'tmax,time,dt,it =',tmax,time,dt,it
      endif
        kdgof = 8
      KVAR=KNODE*KDGOF
        KCOOR=3
        KELEM=32000000
        LR=10000
C        WRITE(*,*) 'KNODE,KDGOF,KVAR,KCOOR,KELEM,LR ='
C        WRITE(*,'(1X,6I7)') KNODE,KDGOF,KVAR,KCOOR,KELEM,LR
      knb1=kdgof*knode*1
      if (knb1/2*2 .lt. knb1) knb1=knb1+1
      kna1=kdgof*knode*2
      kna3=kcoor*knode*2
      kna5=kdgof*knode*2
      kna6=kdgof*knode*2
      kna7=knode*2
      kna8=knode*2
      kna9=knode*2
      kna10=kvar*2
      kna2=kvar*2
      kna4=10000*2
      kna11=kelem*2
      knb2=kelem*1
      if (knb2/2*2 .lt. knb2) knb2=knb2+1
      kna0=1
      kna1=kna1+kna0
      kna2=kna2+kna1
      kna3=kna3+kna2
      kna4=kna4+kna3
      kna5=kna5+kna4
      kna6=kna6+kna5
      kna7=kna7+kna6
      kna8=kna8+kna7
      kna9=kna9+kna8
      kna10=kna10+kna9
      kna11=kna11+kna10
      if (kna11-1.gt.256000000) then
      write(*,*) 'exceed memory of array ia'
      write(*,*) 'memory of ia = 256000000'
      write(*,*) 'memory needed = ',kna11,' in prgram elds3b'
      stop 55555
      endif
      knb0=1
      knb1=knb1+knb0
      knb2=knb2+knb1
      if (knb2-1.gt.128000000) then
      write(*,*) 'exceed memory of array ib'
      write(*,*) 'memory of ib = 128000000'
      write(*,*) 'memory needed = ',knb2,' in prgram elds3b'
      stop 55555
      endif
      call elds3b(knode,kdgof,kvar,kcoor,
     *numtyp,numel,neq,kelem,lr,time,
     *dt,kdeof,ia(kna0),ia(kna1),ia(kna2),
     *ia(kna3),ia(kna4),ia(kna5),ia(kna6),ia(kna7),
     *ia(kna8),ia(kna9),ia(kna10),ib(knb0),ib(knb1),
     *filename)
      return
      end
      subroutine elds3b(knode,kdgof,kvar,kcoor,
     *numtyp,numel,neq,kelem,lr,time,
     *dt,kdeof,u,f,coor,emate,
     *ew,emain,e1a,e2a,e3a,emass,
     *sml,nodvar,node,
     *filename)
      implicit real*8 (a-h,o-z)
      character*40 filename(20)
        DIMENSION NODVAR(KDGOF,KNODE),U(KDGOF,KNODE),COOR(KCOOR,KNODE),
     *ew(kdgof,knode),emain(kdgof,knode),e1a(knode),e2a(knode),
     *e3a(knode),Emass(kvar),
     & F(KVAR),EMATE(10000),SML(KELEM),NODE(KELEM)
 
C ...............................................................C
C         COMPUTE THE SOLUTION BY LEAST SQUARE METHOD            C
C ...............................................................C
C ..... KDGOF NUMBER OF D.O.F
C ..... KNODE NUMBER OF NODES
C ..... NODVAR DENOTE THE EQUATION NUMBER CORRESPONDING THE D.O.F
C ..... U INITIAL VALUE
C ..... COOR COORDINATES
C ..... F RIGHT HAND SIDE TERM OF THE EQUATION
C ..... SML STORE THE ELEMENT STIF,MASS,DAMP MATRIX AND LOAD
C ..... NODE ELEMENT NODAL CONNECTION
C ..... EMATE ELEMENT MATERIAL
C .................................................................
6     FORMAT (1X,26I3)
7     FORMAT (1X,6E12.3)
1001  FORMAT(1X,9I7)
 
C.......OPEN ID0 file
      OPEN (1,FILE=filename(2),FORM='UNFORMATTED',STATUS='unknown')
c        read (1) mid,nid,((NODVAR(I,J),I=1,nid),J=1,KNODE)
      CLOSE (1)
 
C.......COMPUTE NODVAR
c      DO 50 J=1,KNODE
c      DO 30 I=1,KDGOF
c      IF (NODVAR(1,J).lt.-1) then
c      NODVAR(I,J) = nodvar(1,j)
c      else
c      NODVAR(I,J) = 1
c      endif
c30    CONTINUE
c50    CONTINUE
 
      NEQ = 0
      DO 100 J=1,KNODE
      DO 100 I=1,KDGOF
c      IF (NODVAR(I,J).GE.-1) then
      NEQ = NEQ+1
      NODVAR(I,J) = NEQ
c      else
c      N = -NODVAR(I,J)-1
c      NODVAR(I,J) = NODVAR(I,N)
c      endif
      U(I,J) = 0.0
100   CONTINUE
C      WRITE(*,*) 'KDGOF =',KDGOF,' KNODE =',KNODE
C      WRITE (*,*) 'NODVAR ='
C      WRITE (*,6) ((NODVAR(I,J),I=1,KDGOF),J=1,KNODE)
 
C.......OPEN COOR file
      OPEN (1,FILE=filename(3),FORM='UNFORMATTED',STATUS='OLD')
      READ (1) NUMNOD,NCOOR,((COOR(I,J),I=1,NCOOR),J=1,NUMNOD)
      CLOSE(1)
C      WRITE(*,*) 'NUMNOD,NCOOR=',NUMNOD,NCOOR
 
      numtyp = 1
      IF (IT.EQ.0) THEN
      open(11,file='unod',form='unformatted',status='old')
      ELSE
      open(11,file='unod',form='unformatted',status='old')
      ENDIF
      read(11) (e1a(i),i=1,knode),
     &  (e2a(i),i=1,knode),
     &  (e3a(i),i=1,knode)
 
C.......OPEN einform estiff file
      OPEN (1,FILE=filename(4),FORM='UNFORMATTED',STATUS='unknown')
      OPEN (2,FILE=filename(5),FORM='UNFORMATTED',STATUS='unknown')
      DO 110 I=1,NEQ
      Emass(i) = 0.0
110   CONTINUE
      NUMEL=0
C.......OPEN ELEM0 file
      OPEN (3,FILE=filename(6),FORM='UNFORMATTED',STATUS='OLD')
C.......NUMTYP number of element types
      DO 2000 ITYP=1,NUMTYP
C.......READ NODAL CONNECTION POINTS AND MATERIAL No.
C.......NUM .... number of elements
C.......NNODE .... number of nodes for each element
      READ (3) NUM,NNODE,
     *           ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
C      WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODE
C      WRITE(*,*) 'NODE ='
C      WRITE(*,6) ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM)
      NNE = NNODE
      nne = nne-1
      K=0
      DO 116 J=1,NNE
      JNOD = NODE(J)
      IF (JNOD.LE.0) GOTO 116
      DO 115 L=1,KDGOF
      IF (NODVAR(L,JNOD).NE.0) K=K+1
115   CONTINUE
116   CONTINUE
C      WRITE(*,*) 'K =',K
      kk=k*k
      k0=1
      k1=k0+k*k
      k2=k1+k
      k3=k2+k
      k4=k3+k*k
      k5=k4+k*k
      k6=k5+k
      k7=k6+k
      CALL elds3b_1(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE,
     *NUMEL,ITYP,NCOOR,NUM,TIME,DT,NODVAR,COOR,NODE,EMATE,
     &sml(k0),sml(k1),sml(k2),sml(k3),sml(k4),sml(k5),sml(k6),
     &e1a,e2a,e3a,ew,emain,Emass,
     *U)
      NUMEL = NUMEL + NUM
2000  CONTINUE
 
      CLOSE(1)
      CLOSE(2)
      CLOSE(3)
      close(11)
 
C......EMASS LUMPING MATRIX FOR LEAST SQUARE METHOD
      emmax=0.0
      do I=1,NEQ
      if (emmax.lt.emass(I)) emmax=emass(I)
      enddo
      emmin = emmax*1.d-8
      NEQ = 0
      DO 2050 IJ=1,KNODE*KDGOF
      if (emass(IJ).lt.emmin) emass(IJ)=emmin
2050  F(IJ)=0.0D0
      DO 2200 I=1,KNODE
      DO 2100 J=1,KDGOF
      IJ=NODVAR(J,I)
      IF (IJ.LE.0) GOTO 2100
      IF (IJ.GT.NEQ) NEQ = IJ
      F(IJ) = F(IJ)+U(J,I)/EMASS(IJ)
2100  CONTINUE
2200  CONTINUE
      DO 2400 I=1,KNODE
      DO 2300 J=1,KDGOF
      IJ=NODVAR(J,I)
      Ew(J,I) = 0.0
      IF (IJ.LE.0) GOTO 2300
      U(J,I) = F(IJ)
      Ew(J,I) = F(IJ)
2300  CONTINUE
2400  CONTINUE
 
 
C.......OPEN unod FILE AND STORE THE SOLUTION
      OPEN (2,FILE=filename(7),FORM='UNFORMATTED',STATUS='unknown')
      WRITE(2) ((U(J,I),I=1,KNODE),J=1,KDGOF)
      CLOSE (2)
 
C.......OPEN LUMP MATRIX FILE
      OPEN (2,FILE=filename(8),FORM='UNFORMATTED',STATUS='unknown')
      write(2)  (Emass(i),i=1,neq)
      CLOSE (2)
      do 500 inod=1,knode
       call mstress6(6,ew(1,inod),emain(1,inod))
       emain(1,inod)=(emain(4,inod)-emain(6,inod))
       emain(2,inod)=dsqrt(((emain(4,inod)-emain(5,inod))**2
     * +(emain(5,inod)-emain(6,inod))**2
     * +(emain(6,inod)-emain(4,inod))**2)/2.d0)
500   continue
      open(10,file='unods',form='unformatted',status='unknown')
      write(10) (Emain(2,i),i=1,knode),(Emain(4,i),i=1,knode),
     & (Emain(5,i),i=1,knode),(Emain(6,i),i=1,knode)
      close(10)
      END
 
 
      SUBROUTINE elds3b_1(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE,
     & NUMEL,ITYP,NCOOR,NUM,TIME,DT,NODVAR,COOR,NODE,EMATE,
     *es,em,ef,Estifn,Estifv,Emassn,Emassv,
     *e1a,e2a,e3a,ew,emain,Emass,
     & U)
      implicit real*8 (a-h,o-z)
      DIMENSION NODVAR(KDGOF,KNODE),COOR(KCOOR,KNODE),NODE(KELEM),
     & U(KDGOF,KNODE),EMATE(300),
     *es(k,k),em(k),ef(k),e1a(knode),
     *e2a(knode),e3a(knode),ew(kdgof,knode),
     *emain(kdgof,knode),Estifn(k,k),Estifv(kk),Emassn(k),Emassv(k),
     *Emass(1),
     & R(500),PRMT(500),COEF(500),LM(500)
17      FORMAT (1X,15I5)
18      FORMAT (1X,8e9.2)
 
      READ (3) MMATE,NMATE,((EMATE((I-1)*NMATE+J),J=1,NMATE),
     *	I=1,MMATE)
c      WRITE(*,*) 'MMATE =',MMATE,' NMATE =',NMATE
c      WRITE (*,*) 'EMATE ='
c      WRITE (*,18) ((EMATE((I-1)*NMATE+J),J=1,NMATE),
c     *	I=1,MMATE)
        DO 1000 NE=1,NUM
      NR=0
      DO 130 J=1,NNE
      JNOD = NODE((NE-1)*NNODE+J)
        IF (JNOD.LT.0) JNOD = -JNOD
        PRMT(NMATE+7+J) = JNOD
      coef(j+0*nne)=e1a(jnod)
      coef(j+1*nne)=e2a(jnod)
      coef(j+2*nne)=e3a(jnod)
      DO 120 I=1,NCOOR
      NR=NR+1
120   R(NR) = COOR(I,JNOD)
130   CONTINUE
      IMATE = NODE(NNODE*NE)
      DO 140 J=1,NMATE
140   PRMT(J) = EMATE((IMATE-1)*NMATE+J)
      PRMT(NMATE+1)=TIME
      PRMT(NMATE+2)=DT
        PRMT(NMATE+3)=IMATE
        prmt(NMATE+4)=NE
        prmt(NMATE+5)=NUM
        prmt(NMATE+6)=IT
        prmt(NMATE+7)=NMATE
        prmt(NMATE+8)=ITIME
        prmt(NMATE+9)=ITYP
 
      goto 1
1     call bec8g2(r,coef,prmt,es,em,ec,ef,ne)
      goto 2
2     continue
 
C       WRITE(*,*) 'ES EM EF ='
C       DO 555 I=1,K
C555    WRITE(*,18) (ES(I,J),J=1,K)
C       WRITE(*,18) (EM(I),I=1,K)
C       WRITE(*,18) (EF(I),I=1,K)
 
      do 201 i=1,k
      do 201 j=1,k
      Estifn(i,j)=0.0
201   continue
      do 202 i=1,k
      Estifn(i,i)=Estifn(i,i)
      do 202 j=1,k
      Estifn(i,j)=Estifn(i,j)+es(i,j)
202   continue
      do 203 i=1,k
      Emassn(i)=0.0
203   continue
      do 204 i=1,k
      Emassn(i)=Emassn(i)+em(i)
204   continue
 
      L=0
      M=0
      I=0
      DO 700 INOD=1,NNE
      NODI=NODE((NE-1)*NNODE+INOD)
      DO 600 IDGF=1,KDGOF
      INV=NODVAR(IDGF,NODI)
      IF (INV.EQ.0) GOTO 600
      I=I+1
      IF (INV.LT.0) GOTO 305
      L=L+1
      LM(L)=INV
      U(IDGF,NODI)=U(IDGF,NODI)
     *+ef(i)
      Emassv(l)=Emassn(i)
305     J=0
      DO 500 JNOD=1,NNE
      NODJ=NODE((NE-1)*NNODE+JNOD)
      DO 400 JDGF=1,KDGOF
      JNV=NODVAR(JDGF,NODJ)
      IF (JNV.EQ.0) GOTO 400
      J=J+1
 
      IF (JNV.LT.0) GOTO 400
      IF (INV.LT.0) GOTO 310
      M=M+1
      Estifv(m)=Estifn(i,j)
310     CONTINUE
 
 
      IF (INV.LT.0)
     *  U(JDGF,NODJ)=U(JDGF,NODJ)-ESTIFN(I,J)*U(IDGF,NODI)
400     CONTINUE
500     CONTINUE
600     CONTINUE
700     CONTINUE
C       WRITE (*,*) 'U ='
C       WRITE (*,18) ((U(J,I),J=1,KDGOF),I=1,KNODE)
 
      LRD=M
      NER=NUMEL+NE
      DO 800 I=1,L
      J=LM(I)
      Emass(j) = Emass(j) + Emassv(i)
800     CONTINUE
 
1000    CONTINUE
 
      END
